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ABSTRACT 

We present a finite element method (FEM) solver for computation of optical resonance modes in VCSELs. We 
perform a convergence study and demonstrate that high accuracies for 3D setups can be attained on standard 
computers. We also demonstrate simulations of thermo-optical effects in VCSELs. 
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1. INTRODUCTION 

O I Vertical-cavity surface-emitting lasers (VCSELs) are light sources with unique properties and potential applica- 
O . tions of great interest.^ Solution of Maxwell's equations for realistic 3D VCSELs is a challenging task.^ Since 
the VCSEL resonator is realized by distributed Bragg refiectors (DBR), the geometry inherits a pronounced 
multiscale structure. The devices are relatively large (thousands of cubic wavelengths) and contain subwave- 
c/) ■ length DBR layers, very thin active zones and structured apertures. Further, the infinite exterior adjacent to 
the VCSEL containing also a layered structure has to be modeled to obtain realistic predictions of radiation loss 
and lasing threshold. A variety of methods has been used to compute optical VCSEL modes. These include 
FEM, finite difference time-domain methods (FDTD), modal expansion and approximative methods.^^^ In most 
standard approaches the optical problem is restricted to purely ID or to cylindrically symmetric structures. 
Nevertheless, many realistic 3D devices cannot be restricted in this way. This is the reason why reliable full 
3D simulations become so important. Accuracy limitations of state-of-the-art 3D solvers, including also FEM 
solvers, have recently been discussed.^ ^ 

The finite element method is very well suited for simulation of nano-optical systems and devices .l^l^ Its main 
features are the capability of exact geometric modeling due to usage of unstructured meshes and high accuracy at 
low computational cost. The finite element method offers great fiexibility to approximate the solution: different 
mesh sizes h and polynomial ansatz functions of varying degree p can be combined to obtain high convergence 
rates. As a result, very demanding problems can be solved on standard workstations."^ We demonstrate that 
a FEM eigenmode solver with higher-order finite elements, adaptive meshing and a rigorous implementation of 
transparent boundary conditions is a powerful method for 3D VCSEL mode computation. 
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2. MATHEMATICAL BACKGROUND 



The main physical effects in a VCSEL are associated to time scales ranging over several orders of magnitude. 
Since the frequency of the optical modes is much higher than those of all other effects, a time-harmonic ansatz 
for the electric field E{x, z) is well-justified: 

E{x,y,z,t) =e-'^'E{x,y,z), (1) 

where 00 denotes the frequency. Using this ansatz in Maxwell's equations, the following second order equation 
for the electric field can be derived: 

Vx/i-^Vx^ =uj'^eE. (2) 

In this equation no exterior current or charge density sources are present. Physically, the light field of a VCSEL 
is created by coupling of the electron system in the active layer to the eigenmodes of the structure. In Maxwell's 
equations this usually enters via the complex permittivity tensor e (in all relevant optical materials the magnetic 
permeability /i is a constant). The resonance problem then consists of finding pairs (£^,cj), such that Maxwell's 
equations (|2]) on the given geometry are fulfilled. Furthermore, the so called radiation condition has to be 
satisfied which requires that the resonance modes are purely outward radiating. For solving equations (j2j) we 
use the FEM package JCMsuite developed by ZIB and JCMwave.^^ 




(a) (b) 
Figure 1: a) 2D cross section of cylindrically symmetric VCSEL (aperture diameter = 6 ^im, cf., Bienstmann et al.^), the 
inset shows a detail of the finite element triangulation, including part of the 5 nm thin gain region (green), b) Visualization 
of the computed electric field (|^(x, ^, z)|^) of the fundamental mode (A = 980.587 nm, 7 = —2.388 • 10~^). 
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Figure 2: Convergence graphs (reference solution): relative errors of the resonance wavelength A and the damping factor 
7 in dependence on number of unknowns N for different polynomial orders p of the FEM approximation. 



3. VCSEL MODEL AND SIMULATION RESULTS 

The aim of this paper is to demonstrate that very low numerical errors can be reached in full 3D VCSEL 
simulation. In order to quantify the error of a numerical solution we compare it to a reference solution. Because 
for the full 3D VCSEL problem no accurate quantitative results are available as independent reference, we use 
a cylindrically symmetric VCSEL setup in this convergence study. A very accurate solution to this problem can 
be obtained using a 2D solver in cylindrical coordinates.^ Restriction of 3D resonance mode computations to a 
2D cross section due to the cylindrical symmetry leads to substantial savings in computational time and memory 
requirements. The 3D, rotationally symmetric solution can also be compared to results from the literature.'^ 
With the reference solution at hand we perform simulation of the same physical setup using a full 3D FEM 
model. Numerical accuracy of the 3D results from this model is then obtained from comparison to the reference 
solution. 

In a further sub-section, we demonstrate 3D simulation of the thermal lensing effect in a VCSEL. 
3.1 2D reference solution 

As physical model we choose a VCSEL setup as described by Bienstman et ol^ The laser consists of two DBR 
mirrors with alternating AlGaAs-GaAs layers. The InGaAs quantum well layer (gain material) is embedded in 
the central GaAs cavity region. An AlOx aperture is placed in the upper region of the lowest AlGaAs layer 
of the top mirror. The whole structure is situated on a GaAs substrate. Figure [Tal shows a 2D cross section 
through the cylindrically symmetric setup. Note that the layered structure extends infinitely in radial direction 
and modal confinement is reached through the finite AlOx aperture and the finite active region. 

We use a 2D FEM solver in cylindrical coordinates to obtain the near field solution E and the complex 
eigenfrequency u. From u we derive the resonance wavelength A = c • 27r/3?(a;), where c denotes the speed of 
light, and the damping factor 7 = ^(co')/3?(cj), which quantifies the gain/loss of the cavity mode. Figure [Tb] 
shows the electric field intensity distribution of the fundamental VCSEL mode in a pseudo-color representation. 

We have computed solutions to the same physical setup using different numerical resolutions, where we have 
varied finite element degree p and the mesh refinement. With increasing p and increasing mesh refinement also 
the number of unknowns of the FEM problem, A^, increases. Using the solution with highest p and highest 
refinement (p = 6, 4 successive adaptive grid refinements) as reference we can attribute relative numerical errors 
to the computed resonance wavelengths and damping factors. Figure [2] shows how these converge with increasing 
N. We have also checked convergence with respect to the numerical parameters of the perfectly matched layers 



(PML)P We observe that very high accuracies can be reached with relative errors below 10 ^ for the resonance 
wavelength. Note that the relative error of the damping factor is larger due the much smaller total value of the 
imaginary part of the fundamental eigenvalue. However, also here we reach relative accuracies below 0.001%. 

3.2 Full 3D computation 

For full 3D computation of the same VCSEL setup as before we discretize the VCSEL geometry using a 3D 
prismatoidal mesh. Figure [3] (left) shows a visualization of the mesh discretizing the geometry. Note that due 
to mirror symmetries of the device we can restrict the computational domain to one quarter of the full VCSEL, 
on the mirror faces boundary conditions corresponding to the symmetry of the fundamental mode are applied. 
Further, transparent boundary conditions, realized with perfectly matched layers are applied. Figure [3] (right) 
shows a cross section through the computed 3D field distribution {\E{x^ y^z)\'^). As expected the field distribution 
agrees perfectly with the field distribution of the reference solution and the resonance wavelengths agree to high 
precision. 

Again we perform a convergence study where we vary finite element degree p. The relative errors are defined 
as deviations from the 2D cylindrically symmetric reference solution (see above), normalized to the reference 
solution result. Figure H] shows that we reach an accuracy of the resonance wavelength well below 10 ~^ and an 
accuracy of the damping factor below 10~^. Computational times for the full 3D simulations are few minutes 
on a standard workstation (3min {p = 2), 6min {p = 3), resp. 30min {p = 4)). 

In the previous we have simulated a "cold cavity" VCSEL, i.e., without gain in the active zone. Next we 
demonstrate the accuracy of laser threshold computation. For this we perform computations with increasing 
gain in the active zone (imaginary part of the complex refraction index of the quantum well layer, n^). Figure l4bl 
shows the expected linear dependence of the damping factor on material gain n^. We observe that threshold is 
reached at a refractive index Ui ~ 0.011. This value is obtained from 3D simulation at both investigated accuracy 
levels, as well as from the 2D cylindrically symmetric reference simulation. The agreement corresponds to the 
accuracy of the damping factor better than 1% relative deviation as displayed in Fig. HI 

Now, that we have characterized the accuracy level which can be reached by the 3D FEM solver for typ- 
ical VCSEL setups we turn to a "real", i.e., non-rotationally-symmetric VCSEL geometry. Figure [5] shows a 
visualization of the mesh and the fundamental mode of a photonic-crystal VCSEL in a setup as described in.^^ 




Figure 3: 3D layout of the cylindrically symmetric VCSEL with finite element triangulation (left) and visualization of 
the fundamental resonance mode (A = 980.584 nm, 7 = — 2.375 •10-^ right). 
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(a) (b) 
Figure 4: a) Relative error of the resonance wavelength A and the damping factor 7 for 3D simulation versus cylindrical 
2D solution in dependence on number of unknowns N for different polynomial orders p of the FEM approximation, b) 
Threshold gain for 2D and 3D simulations for different polynomial orders p. 




We have used simulation accuracy settings as in the previous results. We therefore expect that the numerical 
accuracy of the wavelength is of the order of 10 ~^ and the numerical accuracy of the damping factor is of the 
order of IQ-^ c.f., Fig.H This is a significantly higher accuracy than reached in the benchmark of Dems et al.^ 
A main difference in our implementation, compared to the FEM implementation in Dems et al.j^is the usage of 
higher-order elements, and of adaptive transparent boundary conditions.'^ 

3.3 Thermo-optical simulations 

Limitations of the output profile quality of VCSELs can be due to unconverted pump power in the active region 
which leads to a local heating and through thermo-optical coupling to inferior laser mode properties. In this 



Figure 6: 3D layout of a VCSEL for thermo-optical simulation with finite element triangulation (left). Visualization of 
the computed temperature profile in an xz-plane (center) and in an xy-plane (right) for AT = Tmax — Tq — 50. 5K in jet 
colormap (Tmax = 350. 5K (dark red), Tmin = 300K (dark blue), Tmax is the maximum temperature in the computational 
domain) . 
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Figure 7: Resonance wavelength of the optical mode. A, in dependence on the maximum temperature difference AT = 
Tmax - To in the VCSEL. 
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Figure 8: a) Electric field intensity along the vertical dimension, b) Electric field intensity along the horizontal dimension, 
subsection we present 3D thermo-optical simulations based on the stationary heat conduction equation 



where k denotes the thermal conductivity, T denotes the temperature and q is the thermal source density. We 
couple this equation to the optical simulation via the thermo-optical correction of the temperature dependent 
relative permittivity e(T) = e(To) + C{T — Tq), where C is the thermo-optical coefficient and Tq is a reference 
temperature. 

We investigate the same VCSEL geometry and material parameters as in the previous Section 13.21 For 
the thermal simulation we need to expand the computational domain (GaAs substrate) in order to set suitable 
boundary conditions, i.e., to fix the boundary temperature of the computational domain, Tq, at 300K. We define 
a heat source in the oxide aperture and in the cavity by setting the thermal source density q in the according 
regions and compute the temperature distribution within the VCSEL by solving the stationary heat equation. 
Fig. [6] shows visualizations of the discretized geometry and of a computed temperature distribution. We then use 
this temperature distribution to define a locally varying permittivity correction throughout the VCSEL geometry, 
and compute the optical cavity modes using the Maxwell solver as in Section 13.21 

We have performed simulations for various heat source powers: Figure [3 shows the almost linear dependence 
of the resonance wavelength A on the maximum temperature difference in the device. Increasing the maximum 
temperature in the device through increasing heat source power causes a refractive index change of the DBR 
structure and a shift of its reflective wavelength. Thus the VCSEL resonance wavelength A is shifted. 

In post-processes, from the near field distribution we exported ID electric field intensity distributions in 
both, vertical and horizontal directions, through the cavity center. The electric field along the vertical dimension 
is nearly independent on temperature variations, as expected, see Fig. [Sal In the horizontal direction we observe 
the constriction of the electric field intensity with increasing temperature, see Fig. [HH Such behaviour is known 
as thermal lensing effect and a limitation for high-power semiconductor lasers. In specific applications, VCSELs 
should be designed for low thermal lensing.^ 



We have presented FEM results for resonance mode computation in VCSELs. Our results for a cylindrically 
symmetric device demonstrate fast convergence and high accuracy of the method. The resonance wavelength and 
the corresponding damping factor can be computed very accurately with relative errors below 10~^, respectively 



V -kVT{x,y,z) = q{x,y,z), 



4. CONCLUSION 



10-^. The reached accuracies are several orders of magnitude higher than the differences between the results 
from different models presented in Ref.*^ 

Comparison of full 3D simulation results to the 2D solution and a detailed convergence analysis confirm the 
high accuracy of the method also for full 3D problems with relative errors of the wavelength and the corresponding 
damping factor below 10~^, respectively 10~^. To our knowledge results on numerical accuracy of 3D VCSEL 
simulations on this level have not previously been published (see, e.g.|^^. We have also presented thermo-optical 
simulations for 3D VCSELs. In future research we plan to investigate more complex 3D device geometries like 
metal cavity surface-emitting microlasers. 
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